###
#
# This is the replication script for "Crisis Management: Personal Financial Well-Being and Public Attitudes Toward Government Intervention"
# Contact:
# Tyler Girard
# tgirard@purdue.edu
#
###

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Set Up HERE ####

library(here)
here()

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Load Packages ####

library(rio)
library(MASS)
library(mice)
library(mitools)
library(tidyverse)
library(modelsummary)
library(table1)
library(ggeffects)
library(ggstance)
library(kableExtra)
library(ggpubr)
library(tictoc)

# for latex output of regression tables
options("modelsummary_format_numeric_latex" = "plain")

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Helper Functions ####

# extract labels
extract_varlabs <- function(data = data){
  #if (!require(expss)) install.packages('expss')
  tmp <- lapply(data, expss::var_lab)
  tmp <- as.data.frame(do.call("rbind", tmp))
  tmp <- tmp %>% rownames_to_column()
  names(tmp) <- c("Variable","Label")
  return(tmp)
}

# Customize modelsummary to extract polr results
tidy_custom.polr <- function(x, ...) {
  s <- lmtest::coeftest(x)
  out <- data.frame(
    term = row.names(s),
    p.value = s[, "Pr(>|t|)"])
  out
}

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Set Seed ####

set.seed(515)

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Load Data ####

# data originally cleaned in Stata
df <- import(here("2 Clean Data", "COVID Survey Clean Data.dta"))


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Extract Labels ####
varlabs <- extract_varlabs(df)



#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Basic Recoding ####

df <- df %>%
  mutate(fed_handling_na = as.factor(fed_handling_na),
         fed_handling = as.factor(fed_handling),
         prov_handling_na = as.factor(prov_handling_na),
         prov_handling = as.factor(prov_handling),
         big = as.factor(big),
         big_na = as.factor(big_na),
         region_na = case_when(Q3_na == 1 ~ "Atlantic",
                               Q3_na == 2 ~ "Quebec",
                               Q3_na == 3 ~ "Ontario",
                               Q3_na == 4 ~ "Prairies",
                               Q3_na == 5 ~ "Alberta",
                               Q3_na == 6 ~ "British Columbia"),
         region_na = fct_relevel(region_na, "Atlantic"),
         region = case_when(Q3 == 1 ~ "Atlantic",
                            Q3 == 2 ~ "Quebec",
                            Q3 == 3 ~ "Ontario",
                            Q3 == 4 ~ "Prairies",
                            Q3 == 5 ~ "Alberta",
                            Q3 == 6 ~ "British Columbia"),
         region = fct_relevel(region, "Atlantic"),
         govt_int = rowSums(select(., bank_relief, rrsp_withdraw, students_money), na.rm = TRUE),
         govt_int = (govt_int/3)-4)


table(df$govt_intervention)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Multiple Imputation ####

df_imp <- df %>%
  select(
    # DVs (need to make govt_int)
    fed_handling_na, prov_handling_na, big_na,
    # IVs
    ppl_struggling_na, fin_hardship_na, cerb, new_debt_na, leftright_na, female, unemployed_na, income_na, children_na, region_na,
    # Index vars
    bank_relief_na, rrsp_withdraw_na, students_money_na,
    # Additional vars for  MI
    Q7_na, #education
    Q8_na #visible minority
  )

# need this dataframe for DV3
df_imp_withmiss <- df_imp %>%
  mutate(govt_int = bank_relief_na + rrsp_withdraw_na + students_money_na,
         govt_int = (govt_int/3)-4)



impvars <- c(
  "fed_handling_na", 
  "prov_handling_na",
  "big_na", 
  "ppl_struggling_na", 
  "fin_hardship_na", 
  "new_debt_na", 
  "leftright_na", 
  "unemployed_na", 
  "income_na", 
  "children_na", 
  "bank_relief_na", 
  "rrsp_withdraw_na", 
  "students_money_na",
  "Q7_na",
  "Q8_na"
)

noimpvars <- c(
  "cerb",
  "female",
  'region_na'
)


names(df_imp)

m0 <-mice(df_imp, maxit = 0)

meth <- m0$method
meth

meth[names(meth) %in% impvars] <- "rf"
meth[names(meth) %in% noimpvars] <- ""

meth

tic()
set.seed(515)
imp <-mice(df_imp, 
           meth = meth,
           print = TRUE, 
           m = 100, 
           maxit = 10,
           seed = 515)
toc()

table(complete(imp)$fed_handling_na)
table(df_imp$fed_handling_na)

table(complete(imp)$leftright_na)
table(df_imp$leftright_na)

#save(imp, file = here("2 Clean Data", "Clean Survey Data (MI Object).RData"))

# to skip re-creating MI data (run time is ~ 10 mins), you can load this file
#load(here("2 Clean Data", "Clean Survey Data (MI Object).RData"))


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ Make index with imputed variables ~ ####

df_imp_mi <- complete(imp, action = "long", include = TRUE)

df_imp_mi <- df_imp_mi %>%
  mutate(
    govt_int = rowSums(select(., bank_relief_na, rrsp_withdraw_na, students_money_na), na.rm = TRUE),
    govt_int = (govt_int/3)-4
  )


df_imp_mi <- as.mids(df_imp_mi)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Descriptive Statistics ####

descriptivetab <- df_imp_withmiss %>%
  mutate(fed_handling_na = case_when(fed_handling_na == 1 ~ "Very bad",
                                     fed_handling_na == 2 ~ "Bad",
                                     fed_handling_na == 3 ~ "Neither",
                                     fed_handling_na == 4 ~ "Good",
                                     fed_handling_na == 5 ~ "Very good"),
         fed_handling_na = fct_relevel(fed_handling_na,
                                       "Very bad", "Bad", "Neither", "Good", "Very good"),
         prov_handling_na = case_when(prov_handling_na == 1 ~ "Very bad",
                                      prov_handling_na == 2 ~ "Bad",
                                      prov_handling_na == 3 ~ "Neither",
                                      prov_handling_na == 4 ~ "Good",
                                      prov_handling_na == 5 ~ "Very good"),
         prov_handling_na = fct_relevel(prov_handling_na,
                                        "Very bad", "Bad", "Neither", "Good", "Very good"),
         big_na = case_when(big_na == 1 ~ "Strongly disagree",
                            big_na == 2 ~ "Disagree",
                            big_na == 3 ~ "Neither",
                            big_na == 4 ~ "Agree",
                            big_na == 5 ~ "Strongly agree"),
         big_na = fct_relevel(big_na,
                              "Strongly disagree", "Disagree", "Neither", "Agree", "Strongly agree"),
         # Predictors
         ppl_struggling_na = case_when(ppl_struggling_na == 1 ~ "Strongly disagree",
                                       ppl_struggling_na == 2 ~ "Disagree",
                                       ppl_struggling_na == 3 ~ "Neither",
                                       ppl_struggling_na == 4 ~ "Agree",
                                       ppl_struggling_na == 5 ~ "Strongly agree"),
         ppl_struggling_na = fct_relevel(ppl_struggling_na,
                                         "Strongly disagree", "Disagree", "Neither", "Agree", "Strongly agree"),
         fin_hardship_na = case_when(fin_hardship_na == 1 ~ "None",
                                     fin_hardship_na == 2 ~ "Little",
                                     fin_hardship_na == 3 ~ "Slight",
                                     fin_hardship_na == 4 ~ "Some",
                                     fin_hardship_na == 5 ~ "Serious"),
         fin_hardship_na = fct_relevel(fin_hardship_na,
                                       "None", "Little", "Slight", "Some", "Serious"),
         cerb = case_when(cerb == 0 ~ "No",
                          cerb == 1 ~ "Yes"),
         cerb = fct_relevel(cerb,
                            "No", "Yes"),
         female = case_when(female == 0 ~ "No",
                            female == 1 ~ "Yes"),
         female = fct_relevel(female,
                              "No", "Yes"),
         unemployed_na = case_when(unemployed_na == 0 ~ "No",
                                   unemployed_na == 1 ~ "Yes"),
         unemployed_na = fct_relevel(unemployed_na,
                                     "No", "Yes"),
         children_na = case_when(children_na == 0 ~ "No",
                                 children_na == 1 ~ "Yes"),
         children_na = fct_relevel(children_na,
                                   "No", "Yes")
  )


# Outcomes
table1::label(descriptivetab$fed_handling_na)<-"Approval of Federal Handling"
table1::label(descriptivetab$prov_handling_na)<-"Approval of Provincial Handling"
table1::label(descriptivetab$govt_int)<-"Approval of Policy Interventions"
table1::label(descriptivetab$big_na)<-"Approval of Basic Income Guarantee"

destab_out <- table1::table1(~fed_handling_na + prov_handling_na + govt_int + big_na, 
                             data=descriptivetab, 
                             topclass="Rtable1-zebra",
                             rowlabelhead="Variable")

table1::t1kable(destab_out, format = "latex")

# Predictors
table1::label(descriptivetab$ppl_struggling_na)<-"Many people are struggling"
table1::label(descriptivetab$fin_hardship_na)<-"I'm feeling financial hardship"
table1::label(descriptivetab$cerb)<-"CERB recipient"
table1::label(descriptivetab$new_debt_na)<-"New debt"
table1::label(descriptivetab$leftright_na)<-"Lef-Right Self-Placement"
table1::label(descriptivetab$female)<-"Female"
table1::label(descriptivetab$unemployed_na)<-"Unemployed/Underemployed"
table1::label(descriptivetab$income_na)<-"Income"
table1::label(descriptivetab$children_na)<-"Children"
table1::label(descriptivetab$region_na)<-"Region"

destab_out2 <- table1::table1(~ppl_struggling_na + fin_hardship_na + 
                                cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na, 
                              data=descriptivetab, 
                              topclass="Rtable1-zebra",
                              rowlabelhead="Variable")

table1::t1kable(destab_out2, format = "latex")



#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#



#### Analysis - Federal Handling ####

# median/mode substitution
m1_fed_mm <- polr(fed_handling ~ ppl_struggling + fin_hardship + 
                    cerb + new_debt + leftright + female + unemployed + income + children + region,
                  data = df,
                  Hess = TRUE)

summary(m1_fed_mm)
nobs(m1_fed_mm)

# listwise deletion
m1_fed_lw <- polr(fed_handling_na ~ ppl_struggling_na + fin_hardship_na + 
                    cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
                  data = df_imp_withmiss,
                  Hess = TRUE)

summary(m1_fed_lw)
nobs(m1_fed_lw)

# multiple imputation
m1_fed_mi <- df_imp_mi %>%
  with(polr(fed_handling_na ~ ppl_struggling_na + fin_hardship_na + 
              cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
            Hess = TRUE
  ))

m1_fed_mi <- pool(m1_fed_mi)
summary(m1_fed_mi)


#### ~ SI File Table ~ ####

m1_results <- list()

m1_results[['Substitution']] <- m1_fed_mm
m1_results[['Listwise Deletion']] <- m1_fed_lw
m1_results[['Multiple Imputation']] <- m1_fed_mi

# Summarize
modelsummary(m1_results,             
             stars = TRUE)


cm <- c( 
  'ppl_struggling' = "Many people are struggling",
  'ppl_struggling_na' = "Many people are struggling",
  "fin_hardship" = "I'm feeling financial hardship",
  "fin_hardship_na" = "I'm feeling financial hardship",
  "cerb" = "CERB",
  "new_debt" = "New Debt",
  "new_debt_na" = "New Debt",
  "leftright" = "Left-Right Self-Placement",
  "leftright_na" = "Left-Right Self-Placement",
  "female" = "Female",
  "unemployed" = "Unemployed/Underemployed",
  "unemployed_na" = "Unemployed/Underemployed",
  "income" = "Income",
  "income_na" = "Income",
  "children" = "Children",
  "children_na" = "Children",
  "regionAlberta" = "Region: Alberta",
  "region_naAlberta" = "Region: Alberta",
  "regionBritish Columbia" = "Region: British Columbia",
  "region_naBritish Columbia" = "Region: British Columbia",
  "regionOntario" = "Region: Ontario",
  "region_naOntario" = "Region: Ontario",
  "regionPrairies" = "Region: Prairies",
  "region_naPrairies" = "Region: Prairies",
  "regionQuebec" = "Region: Quebec",
  "region_naQuebec" = "Region: Quebec"
)


tab_m1 <- modelsummary(m1_results,
                       title = 'Approval of Federal Handling',
                       coef_map = cm,
                       estimate = "{estimate} ({std.error}){stars}",
                       statistic = NULL,
                       gof_omit = "edf")

#save_kable(tab_m1, file = here("4 Tables", "SI_RegTab_Federal.png"))

modelsummary(m1_results,
             title = 'Approval of Federal Handling',
             coef_map = cm,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             gof_omit = "edf",
             output = "latex")


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#



#### Analysis - Prov Handling ####

# median/mode substitution
m2_prov_mm <- polr(prov_handling ~ ppl_struggling + fin_hardship + 
                     cerb + new_debt + leftright + female + unemployed + income + children + region,
                   data = df,
                   Hess = TRUE)

summary(m2_prov_mm)
nobs(m2_prov_mm)

# listwise deletion
m2_prov_lw <- polr(prov_handling_na ~ ppl_struggling_na + fin_hardship_na + 
                     cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
                   data = df_imp_withmiss,
                   Hess = TRUE)

summary(m2_prov_lw)
nobs(m2_prov_lw)

# multiple imputation
m2_prov_mi <- df_imp_mi %>%
  with(polr(prov_handling_na ~ ppl_struggling_na + fin_hardship_na + 
              cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
            Hess = TRUE
  ))

m2_prov_mi <- pool(m2_prov_mi)
summary(m2_prov_mi)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### ~ SI File Table ~ ####

m2_results <- list()

m2_results[['Substitution']] <- m2_prov_mm
m2_results[['Listwise Deletion']] <- m2_prov_lw
m2_results[['Multiple Imputation']] <- m2_prov_mi

# Summarize
modelsummary(m2_results,             
             stars = TRUE)


cm <- c( 
  'ppl_struggling' = "Many people are struggling",
  'ppl_struggling_na' = "Many people are struggling",
  "fin_hardship" = "I'm feeling financial hardship",
  "fin_hardship_na" = "I'm feeling financial hardship",
  "cerb" = "CERB",
  "new_debt" = "New Debt",
  "new_debt_na" = "New Debt",
  "leftright" = "Left-Right Self-Placement",
  "leftright_na" = "Left-Right Self-Placement",
  "female" = "Female",
  "unemployed" = "Unemployed/Underemployed",
  "unemployed_na" = "Unemployed/Underemployed",
  "income" = "Income",
  "income_na" = "Income",
  "children" = "Children",
  "children_na" = "Children",
  "regionAlberta" = "Region: Alberta",
  "region_naAlberta" = "Region: Alberta",
  "regionBritish Columbia" = "Region: British Columbia",
  "region_naBritish Columbia" = "Region: British Columbia",
  "regionOntario" = "Region: Ontario",
  "region_naOntario" = "Region: Ontario",
  "regionPrairies" = "Region: Prairies",
  "region_naPrairies" = "Region: Prairies",
  "regionQuebec" = "Region: Quebec",
  "region_naQuebec" = "Region: Quebec"
)


tab_m2 <- modelsummary(m2_results,
                       title = 'Approval of Provincial Handling',
                       coef_map = cm,
                       estimate = "{estimate} ({std.error}){stars}",
                       statistic = NULL,
                       gof_omit = "edf")

#save_kable(tab_m2, file = here("4 Tables", "SI_RegTab_Provincial.png"))

modelsummary(m2_results,
             title = 'Approval of Provincial Handling',
             coef_map = cm,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             gof_omit = "edf",
             output = "latex")



#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### Analysis - Policy Interventions ####

# median/mode substitution
m3_policy_mm <- lm(govt_int ~ ppl_struggling + fin_hardship + 
                     cerb + new_debt + leftright + female + unemployed + income + children + region,
                   data = df)

summary(m3_policy_mm)
nobs(m3_policy_mm)

# listwise deletion
m3_policy_lw <- lm(govt_int ~ ppl_struggling_na + fin_hardship_na + 
                     cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
                   data = df_imp_withmiss)

summary(m3_policy_lw)
nobs(m3_policy_lw)

# multiple imputation
m3_policy_mi <- df_imp_mi %>%
  with(lm(govt_int ~ ppl_struggling_na + fin_hardship_na + 
            cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na
  ))

m3_policy_mi <- pool(m3_policy_mi)
summary(m3_policy_mi)


#### ~ SI File Table ~ ####

m3_results <- list()

m3_results[['Substitution']] <- m3_policy_mm
m3_results[['Listwise Deletion']] <- m3_policy_lw
m3_results[['Multiple Imputation']] <- m3_policy_mi

# Summarize
modelsummary(m3_results,             
             stars = TRUE)


cm <- c( 
  'ppl_struggling' = "Many people are struggling",
  'ppl_struggling_na' = "Many people are struggling",
  "fin_hardship" = "I'm feeling financial hardship",
  "fin_hardship_na" = "I'm feeling financial hardship",
  "cerb" = "CERB",
  "new_debt" = "New Debt",
  "new_debt_na" = "New Debt",
  "leftright" = "Left-Right Self-Placement",
  "leftright_na" = "Left-Right Self-Placement",
  "female" = "Female",
  "unemployed" = "Unemployed/Underemployed",
  "unemployed_na" = "Unemployed/Underemployed",
  "income" = "Income",
  "income_na" = "Income",
  "children" = "Children",
  "children_na" = "Children",
  "regionAlberta" = "Region: Alberta",
  "region_naAlberta" = "Region: Alberta",
  "regionBritish Columbia" = "Region: British Columbia",
  "region_naBritish Columbia" = "Region: British Columbia",
  "regionOntario" = "Region: Ontario",
  "region_naOntario" = "Region: Ontario",
  "regionPrairies" = "Region: Prairies",
  "region_naPrairies" = "Region: Prairies",
  "regionQuebec" = "Region: Quebec",
  "region_naQuebec" = "Region: Quebec",
  "(Intercept)" = "Intercept"
)


tab_m3 <- modelsummary(m3_results,
                       title = 'Approval of Policy Interventions',
                       coef_map = cm,
                       estimate = "{estimate} ({std.error}){stars}",
                       statistic = NULL,
                       gof_omit = "edf")

#save_kable(tab_m3, file = here("4 Tables", "SI_RegTab_Policy.png"))

modelsummary(m3_results,
             title = 'Approval of Policy Interventions',
             coef_map = cm,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             gof_omit = "edf",
             output = "latex")


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Analysis - Basic Income ####


# median/mode substitution
m4_big_mm <- polr(big ~ ppl_struggling + fin_hardship + 
                    cerb + new_debt + leftright + female + unemployed + income + children + region,
                  data = df,
                  Hess = TRUE)

summary(m4_big_mm)
nobs(m4_big_mm)

# listwise deletion
m4_big_lw <- polr(big_na ~ ppl_struggling_na + fin_hardship_na + 
                    cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
                  data = df_imp_withmiss,
                  Hess = TRUE)

summary(m4_big_lw)
nobs(m4_big_lw)

# multiple imputation
m4_big_mi <- df_imp_mi %>%
  with(polr(big_na ~ ppl_struggling_na + fin_hardship_na + 
              cerb + new_debt_na + leftright_na + female + unemployed_na + income_na + children_na + region_na,
            Hess = TRUE
  ))

m4_big_mi <- pool(m4_big_mi)
summary(m4_big_mi)

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ SI File Table ~ ####

m4_results <- list()

m4_results[['Substitution']] <- m4_big_mm
m4_results[['Listwise Deletion']] <- m4_big_lw
m4_results[['Multiple Imputation']] <- m4_big_mi

# Summarize
modelsummary(m4_results,             
             stars = TRUE)

cm <- c( 
  'ppl_struggling' = "Many people are struggling",
  'ppl_struggling_na' = "Many people are struggling",
  "fin_hardship" = "I'm feeling financial hardship",
  "fin_hardship_na" = "I'm feeling financial hardship",
  "cerb" = "CERB",
  "new_debt" = "New Debt",
  "new_debt_na" = "New Debt",
  "leftright" = "Left-Right Self-Placement",
  "leftright_na" = "Left-Right Self-Placement",
  "female" = "Female",
  "unemployed" = "Unemployed/Underemployed",
  "unemployed_na" = "Unemployed/Underemployed",
  "income" = "Income",
  "income_na" = "Income",
  "children" = "Children",
  "children_na" = "Children",
  "regionAlberta" = "Region: Alberta",
  "region_naAlberta" = "Region: Alberta",
  "regionBritish Columbia" = "Region: British Columbia",
  "region_naBritish Columbia" = "Region: British Columbia",
  "regionOntario" = "Region: Ontario",
  "region_naOntario" = "Region: Ontario",
  "regionPrairies" = "Region: Prairies",
  "region_naPrairies" = "Region: Prairies",
  "regionQuebec" = "Region: Quebec",
  "region_naQuebec" = "Region: Quebec"
)


tab_m4 <- modelsummary(m4_results,
                       title = 'Approval of Basic Income Guarantee',
                       coef_map = cm,
                       estimate = "{estimate} ({std.error}){stars}",
                       statistic = NULL,
                       gof_omit = "edf")

#save_kable(tab_m4, file = here("4 Tables", "SI_RegTab_BIG.png"))

modelsummary(m4_results,
             title = 'Approval of Basic Income Guarantee',
             coef_map = cm,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             gof_omit = "edf",
             output = "latex")



#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### Figures ####

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ Coefficient Plots ~ ####

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ ~ Govt Approval ~ ~ ####


cplot_plotdat1 <- broom::tidy(m1_fed_lw, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate(model = "Federal Approval")

cplot_plotdat2 <- broom::tidy(m2_prov_lw, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate(model = "Provincial Approval")



cplot_plotdat <- rbind(cplot_plotdat1,cplot_plotdat2) %>%
  mutate(term = case_when(term == "ppl_struggling" | term == "ppl_struggling_na" ~ "Many people are struggling",
                          term == "fin_hardship" | term == "fin_hardship_na" ~ "I'm feeling financial hardship",
                          term == "cerb" ~ "CERB",
                          term == "new_debt" | term == "new_debt_na" ~ "New Debt",
                          term == "leftright" | term == "leftright_na" ~ "Left-Right Self-Placement",
                          term == "female" ~ "Female",
                          term == "unemployed" | term == "unemployed_na" ~ "Unemployed/Underemployed",
                          term == "income" | term == "income_na" ~ "Income",
                          term == "children" | term == "children_na" ~ "Children"),
         newterm = fct_relevel(term,
                               "Children",
                               "Income",
                               "Unemployed/Underemployed",
                               "Female",
                               "Left-Right Self-Placement",
                               "New Debt",
                               "CERB",
                               "I'm feeling financial hardship",
                               "Many people are struggling"),
         plotcol = case_when(conf.low >= 0 | conf.high <=0 ~ "(95% CI Does Not Contain 0)",
                             TRUE ~ "(95% CI Contains 0)"),
         plotcolcomb = paste(model, plotcol)) %>%
  drop_na(term)

coefplot <- ggplot(cplot_plotdat, aes(shape = model, x = estimate, y = newterm, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0, col = "indianred",lty=4) +
  geom_pointrangeh(aes(linetype = model), position = position_dodge(width = 0.5)) +
  scale_linetype_manual(values = c(5,1)) +
  scale_shape_manual(values = c(5,19)) +
  scale_x_continuous(limits = c(-0.9,0.9)) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(x = NULL, 
       y = NULL,
       title = NULL)

coefplot


#ggsave(coefplot, filename = here("3 Figures", "Manuscript_Govt Approval Coef Plot.png"), scale = 0.5, dpi = 600)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ ~ Policy + BIG ~ ~ ####

cplot_plotdat3 <- broom::tidy(m3_policy_lw, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate(model = "Policy Intervention")

cplot_plotdat4 <- broom::tidy(m4_big_lw, conf.int = TRUE) %>%
  select(term, estimate, conf.low, conf.high) %>%
  mutate(model = "Basic Income Guarantee")



cplot_plotdat2 <- rbind(cplot_plotdat3,cplot_plotdat4) %>%
  mutate(term = case_when(term == "ppl_struggling" | term == "ppl_struggling_na" ~ "Many people are struggling",
                          term == "fin_hardship" | term == "fin_hardship_na" ~ "I'm feeling financial hardship",
                          term == "cerb" ~ "CERB",
                          term == "new_debt" | term == "new_debt_na" ~ "New Debt",
                          term == "leftright" | term == "leftright_na" ~ "Left-Right Self-Placement",
                          term == "female" ~ "Female",
                          term == "unemployed" | term == "unemployed_na" ~ "Unemployed/Underemployed",
                          term == "income" | term == "income_na" ~ "Income",
                          term == "children" | term == "children_na" ~ "Children"),
         newterm = fct_relevel(term,
                               "Children",
                               "Income",
                               "Unemployed/Underemployed",
                               "Female",
                               "Left-Right Self-Placement",
                               "New Debt",
                               "CERB",
                               "I'm feeling financial hardship",
                               "Many people are struggling"),
         plotcol = case_when(conf.low >= 0 | conf.high <=0 ~ "(95% CI Does Not Contain 0)",
                             TRUE ~ "(95% CI Contains 0)"),
         plotcolcomb = paste(model, plotcol)) %>%
  drop_na(term)

coefplot2 <- ggplot(cplot_plotdat2, aes(shape = model, x = estimate, y = newterm, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0, col = "indianred",lty=4) +
  geom_pointrangeh(aes(linetype = model), position = position_dodge(width = 0.5)) +
  scale_linetype_manual(values = c(5,1)) +
  scale_shape_manual(values = c(5,19)) +
  scale_x_continuous(limits = c(-0.9,0.9)) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(x = NULL, 
       y = NULL,
       title = NULL)

coefplot2


#ggsave(coefplot2, filename = here("3 Figures", "Manuscript_Policy + BIG Coef Plot.png"), scale = 0.5, dpi = 600)

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### ~ Predicted Probs Plots ~ ####

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### ~ ~ Fed Approval ~ ~ ####

ppdat1_1 <- ggeffect(m1_fed_lw, terms = ("ppl_struggling_na")) %>%
  mutate(model = "Many people are struggling",
         response.level = case_when(response.level == "X1" ~ "Very Bad",
                                    response.level == "X2" ~ "Bad",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Good",
                                    response.level == "X5" ~ "Very Good"),
         fct_relevel(response.level,
                     "Very Bad", "Bad", "Neither", "Good", "Very Good"),
         plotord = case_when(response.level == "Very Bad" ~ 1,
                             response.level == "Bad" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Good" ~ 4,
                             response.level == "Very Good" ~ 5))


ppdat1_2 <- ggeffect(m1_fed_lw, terms = ("fin_hardship_na")) %>%
  mutate(model = "I'm feeling financial hardship",
         response.level = case_when(response.level == "X1" ~ "Very Bad",
                                    response.level == "X2" ~ "Bad",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Good",
                                    response.level == "X5" ~ "Very Good"),
         fct_relevel(response.level,
                     "Very Bad", "Bad", "Neither", "Good", "Very Good"),
         plotord = case_when(response.level == "Very Bad" ~ 1,
                             response.level == "Bad" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Good" ~ 4,
                             response.level == "Very Good" ~ 5))


ppdat1 <- rbind(ppdat1_1, ppdat1_2)

pp_m1 <- ggplot(ppdat1, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = model)) +
  geom_ribbon(aes(fill = model), alpha = 0.2) +
  scale_color_manual(values = c("#D16103", "#4E84C4")) +
  scale_fill_manual(values = c("#D16103", "#4E84C4")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~fct_reorder(response.level, plotord), ncol = 3) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(
    x = NULL,
    y = "Predicted Probability"
  )

pp_m1

#ggsave(pp_m1, filename = here("3 Figures", "Manuscript_Predicted Probs for Fed Approval.png"))



#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### ~ ~ Prov Approval ~ ~ ####

ppdat2_1 <- ggeffect(m2_prov_lw, terms = ("ppl_struggling_na")) %>%
  mutate(model = "Many people are struggling",
         response.level = case_when(response.level == "X1" ~ "Very Bad",
                                    response.level == "X2" ~ "Bad",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Good",
                                    response.level == "X5" ~ "Very Good"),
         fct_relevel(response.level,
                     "Very Bad", "Bad", "Neither", "Good", "Very Good"),
         plotord = case_when(response.level == "Very Bad" ~ 1,
                             response.level == "Bad" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Good" ~ 4,
                             response.level == "Very Good" ~ 5))


ppdat2_2 <- ggeffect(m2_prov_lw, terms = ("fin_hardship_na")) %>%
  mutate(model = "I'm feeling financial hardship",
         response.level = case_when(response.level == "X1" ~ "Very Bad",
                                    response.level == "X2" ~ "Bad",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Good",
                                    response.level == "X5" ~ "Very Good"),
         fct_relevel(response.level,
                     "Very Bad", "Bad", "Neither", "Good", "Very Good"),
         plotord = case_when(response.level == "Very Bad" ~ 1,
                             response.level == "Bad" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Good" ~ 4,
                             response.level == "Very Good" ~ 5))


ppdat2 <- rbind(ppdat2_1, ppdat2_2)

pp_m2 <- ggplot(ppdat2, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = model)) +
  geom_ribbon(aes(fill = model), alpha = 0.2) +
  scale_color_manual(values = c("#D16103", "#4E84C4")) +
  scale_fill_manual(values = c("#D16103", "#4E84C4")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~fct_reorder(response.level, plotord), ncol = 3) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(
    x = NULL,
    y = "Predicted Probability"
  )

pp_m2

#ggsave(pp_m2, filename = here("3 Figures", "Manuscript_Predicted Probs for Prov Approval.png"))


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### ~ ~ ~ Combine Fed and Prov ~ ~ ~ ####

pp_m1_comb <- ggplot(ppdat1, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = model)) +
  geom_ribbon(aes(fill = model), alpha = 0.2) +
  scale_color_manual(values = c("#D16103", "#4E84C4")) +
  scale_fill_manual(values = c("#D16103", "#4E84C4")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~fct_reorder(response.level, plotord), ncol = 5) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(
    x = NULL,
    y = "Predicted Probability",
    title = "Approval of Federal Handling of Crisis"
  )

pp_m2_comb <- ggplot(ppdat2, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = model)) +
  geom_ribbon(aes(fill = model), alpha = 0.2) +
  scale_color_manual(values = c("#D16103", "#4E84C4")) +
  scale_fill_manual(values = c("#D16103", "#4E84C4")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~fct_reorder(response.level, plotord), ncol = 5) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(
    x = NULL,
    y = "Predicted Probability",
    title = "Approval of Provincial Handling of Crisis"
  )


pp_approval_comb <- ggarrange(pp_m1_comb, pp_m2_comb, ncol = 1,
                              #labels = list("Approval of Federal Handling of Crisis", "Approval of Provincial Handling of Crisis"),
                              legend = "bottom",
                              common.legend = TRUE)

#ggsave(pp_approval_comb, filename = here("3 Figures", "Manuscript_Predicted Probs for Fed + Prov Approval.png"), dpi = 600)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#


#### ~ ~ Basic Income ~ ~ ####

ppdat4_1 <- ggeffect(m4_big_lw, terms = ("ppl_struggling_na")) %>%
  mutate(model = "Many people are struggling",
         response.level = case_when(response.level == "X1" ~ "Strongly Disagree",
                                    response.level == "X2" ~ "Disagree",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Agree",
                                    response.level == "X5" ~ "Strongly Agree"),
         fct_relevel(response.level,
                     "Strongly Disagree", "Disagree", "Neither", "Agree", "Strongly Agree"),
         plotord = case_when(response.level == "Strongly Disagree" ~ 1,
                             response.level == "Disagree" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Agree" ~ 4,
                             response.level == "Strongly Agree" ~ 5))


ppdat4_2 <- ggeffect(m4_big_lw, terms = ("fin_hardship_na")) %>%
  mutate(model = "I'm feeling financial hardship",
         response.level = case_when(response.level == "X1" ~ "Strongly Disagree",
                                    response.level == "X2" ~ "Disagree",
                                    response.level == "X3" ~ "Neither",
                                    response.level == "X4" ~ "Agree",
                                    response.level == "X5" ~ "Strongly Agree"),
         fct_relevel(response.level,
                     "Strongly Disagree", "Disagree", "Neither", "Agree", "Strongly Agree"),
         plotord = case_when(response.level == "Strongly Disagree" ~ 1,
                             response.level == "Disagree" ~ 2,
                             response.level == "Neither" ~ 3,
                             response.level == "Agree" ~ 4,
                             response.level == "Strongly Agree" ~ 5))

ppdat4 <- rbind(ppdat4_1, ppdat4_2)

pp_m4 <- ggplot(ppdat4, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = model)) +
  geom_ribbon(aes(fill = model), alpha = 0.2) +
  scale_color_manual(values = c("#D16103", "#4E84C4")) +
  scale_fill_manual(values = c("#D16103", "#4E84C4")) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~fct_reorder(response.level, plotord), ncol = 3) +
  theme_bw() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank()
  ) +
  labs(
    x = NULL,
    y = "Predicted Probability"
  )

pp_m4

#ggsave(pp_m4, filename = here("3 Figures", "Manuscript_Predicted Probs for Basic Income.png"), dpi = 600)


#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

#### End of Script ####

#---------------------------------------------------------------------------------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------------------------------------------------------------------------------#

